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INTRODUCTION 


This report summarizes the progress made under the NASA Grant NAG-3- 
768 titled "Numerical study of the effects of icing on viscous flow over wings". 
The work was carried out by Dr. L. N. Sankar, the principal investigator, and Dr. 
Oh J. Kwon with the assistance of the following graduate research assistants: Mr. 
Jiunn-Chi Wu, Mr. Ashok Bangalore and Mr. Napporn Phaengsook. Another 
student, Mr. Olympio Melio, not supported under this project, also contributed to 
the work reported. 


The research effort lead to the development of 2-D and 3-D computational 
tools for the prediction of viscous flow over iced wings and airfoils. Much of the 
work has already been published. Here is the list of all the publications, 
supported by the present grant. 

Refereed Publications: 

18. Wu, Jiunn-Chi, Huff, D. and Sankar, L.N., "Evaluation of Three 
Turbulence Models in Static Airloads and Dynamic Stall Predictions," 
Journal of Aircraft . Vol. 27, No. 4, April 1990, pp382-384. 


Invited Conference Keynote Presentation 

1. Sankar, L. N., Kwon, O. J., Bangalore, A., Phaengsook, N. and 
Mello, O., "Effects of Icing on the Performance of Lifting Surfaces," Invited 
Lecture, Workshop on Aircraft Icing and Transition, Ecole Polytechnique, 
University of Montreal, Montreal, Canada, September 20-21, 1993. 

Conference Presentations with Proceedings (Non-Refereed) 


1. Potapczuk, M. G., Bragg, M. B. , Kwon, O. J. and Sankar, L. N., 
"Simulation of Iced Wing Aerodynamics," Proceedings of the AGARD 
Conference on Effects of Adverse Weather on Aerodynamics, AGARD 
CP-496, April 29 - May 1 , 1 991 . 

Conference Presentations without Proceedings (Non-Refereed) 

1. Wu, Jiunn-Chi, Huff, D. and Sankar, L. N., "A Comparison of Three 
Turbulence Models for the Prediction of Steady and Unsteady Airloads," 
AIAA Paper 89-0609. 



2. Kwon, O. J. and Sankar, L. N., "Numerical Simulation of the Effects 
of Icing on the Aerodynamic Characteristics of a Rectangular Wing," AIAA 
28th Aerospace Sciences Meeting, Reno, Nevada, January 1990. 


3. Kwon, O. J. and Sankar, L. N., "Numerical Study of the Effects of 
Icing on the Hover Performance of Rotorcraft," AIAA Paper 91-0662, 
January 1991. 

4. Kwon, O. J. and Sankar, L. N., "Numerical Investigation of 
Performance Degradation of Wings and Rotors due to icing," AIAA Paper 
92-0412. 

5. Sankar, L. N., Phaengsook, N. and Bangalore, A., "Effects of icing 
on the Aerodynamic Performance of High Lift Airfoils," AIAA Paper 93- 
0026. 

6. Mello, O. A. F. and Sankar, L. N., "A Hybrid Navier-Stokes/Full 
Potential Method for the Prediction of Iced Wing Aerodynamics," AIAA 
Paper 94-0489 

7. Bangalore, A., Phaengsook, N. and Sankar, L. N., "Application of a 
Third Order Upwind Scheme to Viscous Flow over Clean and Iced Winqs " 
AIAA Paper 94-0485. 


During the final year of this grant ((January 1 , 1 994 - December 31 , 
1994), work was completed on an efficient hybrid procedure that may be used to 
study clean and iced wing aerodynamics. This work resulted in the Ph. D. 
dissertation of Mr. Olympio Mello. A draft copy of Mr. Mello's Ph. D. thesis 
dissertation is enclosed as an appendix. A revised copy of Mr. Mello's 
dissertation will be mailed to the sponsor in January 1995, after Mr. Mello 
successfully defends his thesis work. 
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SUMMARY 


An improved hybrid method for computing unsteady compressible viscous flows 
is presented. This method divides the computational domain into two zones. In the 
outer zone, the unsteady full-potential equation (FPE) is solved. In the inner zone, the 
Navier-Stokes equations are solved using a diagonal form of an alternating-direction 
implicit (ADI) approximate factorization procedure. The two zones are tightly coupled 
so that steady and unsteady flows may be efficiently solved. Characteristic-based 
viscous/inviscid interface boundary conditions are employed to avoid spurious 
reflections at that interface. The resulting CPU times are less than 60% of the required 
for a full-blown Navier-Stokes analysis for steady flow applications and about 60% of 
the Navier-Stokes CPU times for unsteady flows in non-vector processing machines. 
Applications of the method are presented for a rectangular NACA 0012 wing in low 
subsonic steady flow at moderate and high angles of attack, and for a F-5 wing in 
steady and unsteady subsonic and transonic flows. Steady surface pressures are in very 
good agreement with experimental data and are essentially identical to Navier-Stokes 
predictions. Density contours show that shocks cross the viscous/inviscid interface 
smoothly, so that the accuracy of full Navier-Stokes equations can be retained with a 
significant savings in computational time. 
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CHAPTER I 


INTRODUCTION 


1.1. Overview 

The development of computational fluid dynamics has brought to the industry 
and research communities a variety of methodologies based on the transonic small 
disturbance equation (TSD), full-potential equation (FPE), Euler equations and Navier- 
Stokes equations 1 . TSD- and FPE-based methods have been extensively used to 
compute complex configurations. These methods, in some cases, have been coupled to 
interactive boundary layer analyses to allow solution of problems where viscous effects 
can be included in a limited way. 

For problems where substantial separation occurs, the TSD and FPE techniques 
coupled with interactive boundary-layer analysis are not adequate, since the concept of 
a boundary layer is no longer applicable. For these cases, Navier-Stokes methods are 
clearly needed. However, these are still computationally expensive, and have seen 
limited practical use for complete configurations due to this factor. This becomes 
especially evident for problems where extensive computations are needed, such as the 
prediction of transonic flutter 2 . 
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Hie present method is an extension of the work initiated by Sankar et al. 1 , who 
developed a zonal Navier-Stokes/Full-Potential solver, which was subsequently 
extended to rotors by Tsung and Sankar^ The approach used here is to solve the Full- 
Potential Equation in an outer zone, away from solid surfaces and viscous regions, and 
solve the Navier-Stokes equations in an inner zone, where viscous effects are essential 
This approach is schematically illustrated in Fig. 1.1. This results in a highly efficient 
solver that retains the accuracy of the Navier-Stokes methodology near the solid 
surface, and the simplicity of a potential flow solver away from solid surfaces. 



Fig. 1.1: Partitioning of Computational Domain into Inner and Outer Zones. 
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The Full-Potential solver used in the outer zone solves the unsteady 
compressible FPE in strong conservation form using the artificial compressibility 
concept and employing a strongly implicit procedure 4 * 5 . The Navier-Stokes solver used 
in the inner zone was developed by Sankar et al. and extensively tested in a variety of 
problems 6 " 12 . It employs a diagonal form 13 of an alternating-direction implicit (ADI) 
approximate factorization procedure 14 . 

Historically, coupling potential flow to viscous flow via boundary-layer 
analysis has proved troublesome at the separation point and in the recirculation region. 
Since we are computing the full Navier-Stokes equations in time-dependent form in the 
inner region, the above difficulties associated with boundary-layer methods are 
avoided. 


1.2. Unsteady Transonic Flow an d Aeroelastic Problems 

Transonic flow is characterized by the presence of regions of supersonic flow 
embedded in a subsonic region, as illustrated in Fig. 1.2. Mathematically, the governing 
equations are inherently non-linear, a fact that has prevented the application of 
traditional analytical tools and early numerical methods to the analysis of such a flow 
condition. In addition, transonic flows tend to be more unsteady and three-dimensional 
than purely subsonic and supersonic flows 15 . Despite these difficulties, flight in the 
transonic range is highly desirable for commercial airplanes which achieve their best 
cruise performance at transonic speeds 16 . This flow regime is also encountered by 
modem high performance aircraft during maneuvers 17 , helicopter rotor blades 18 , 
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turbomachinery 19 , launch vehicles in their initial stages of flight 20 re-entry bodies at 
hypersonic speeds and even bluff bodies at subsonic speeds 21 . 



Fig. 1.2: Mixed Flow Regions in Transonic Flow. 

In non-steady flow situations. The presence of a supersonic region embedded in 
a subsonic region causes downstream disturbances to be propagated upstream with a 
considerable time lag, which results in significant out-of-phase forces. 

It has been known for quite some time 22 > 2 ^ that transonic flow conditions are 
critical for flutter, with the flutter dynamic pressure being substantially reduced for 
Mach numbers near unity, in a phenomenon that has been called “transonic dip” 2 ’ 24 . 
This problem is illustrated in Fig. 1.3, from Ref. 2. The severity of flutter at transonic 
speeds is linked to the presence of moving shock waves over the wing surface 2 ^. 
Tijdeman 24 * 2 ^ identified the following types of shock motion: 
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Flutter Dynamic Pressure 


Sinusoidal Shock- Wave Motion (Type A): The shock moves almost 
sinusoidally and remains present during the complete cycle of oscillation 
although its strength varies. Due to the dynamic effect, phase shifts exist 
between the model motion and shock position and between shock strength 
and shock position. The maximum shock strength is not reached during the 
maximum downstream position of the shock, as in quasi-steady flow, but 
during its upstream motion. 

Interrupted Shock-Wave Motion (Type B): This motion is similar to Type A, 
but now the magnitude of the periodic change in shock strength becomes 
larger than the mean steady shock strength and, as a consequence, the shock 
wave disappears during a part of its backward motion. 

Upstream-Propagated Shock Waves (Type C): At slightly supercritical Mach 
number a third type of periodic shock-wave motion is observed, which 
differs completely from the preceding types. Periodically a shock wave is 
formed on die upper surface of the airfoil. This shock moves upstream while 
increasing its strength, The shock wave weakens again, but continues its 
upstream motion, leaves the airfoil from the leading edge, and propagates 
upstream into the oncoming flow as a (weak) free shock wave. This 

phenomenon is repeated periodically and alternates between upper and lower 
surface.” 



Fig. 1.3: Flutter Dynamic Pressure Variation with Mach Number (after Ref. 2). 
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From the above considerations, it is clear that accurate flutter predictions 
depend on the ability of the computational fluid dynamics procedure to predict correct 
shock strength and location, in a time-accurate fashion. Other aeroelastic problems, 
such as tail buffet 2 ’ 17 , are more demanding and require advanced turbulence models, 
since significant separation is characteristic to this phenomenon. 


1.3. Historical Perspective 

Numerical computation of unsteady transonic flows has been one of the major 
challenges to aerodynamicists. Consequently, progress in this area has followed closely 
the development of Computational Fluid Dynamics (CFD). Reviews of the state of the 
art have been published every few years 15,27-30,17. AGARD conferences have been the 
focus of much of the pertinent work, and they have had specific reviews 31-33 . 

Early studies on unsteady transonic flow were impractical because of the 
inherent nonlinearity of that flow range, which prevented the use of the available 
analytical tools, such as panel and doublet- and vortex-lattice methods35>36^ § 0 
transonic flutter predictions had to rely on experiments 22 ’ 34 . 

The numerical computation of transonic flows was initiated by the pioneering 
work of Murman and Cole 37 , which gave rise to substantial development in methods 
for solution of the Transonic Small Disturbance (TSD) equation. Early attempts at 
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numerical solution on harmonically oscillating airfoils 38 were obtained in a two-step 
process, first solving the steady nonlinear problem using the steady transonic small 
disturbance equation and then solving a linear perturbation equation for small 
perturbations about this steady solution. 

Time-accurate solutions to the TSD equation were obtained by Ballhaus et 
al.39,40 xhi s formulation was used by Ballhaus and Gooijian 41 to construct the code 
LTRAN2, which would later be subsequently improved and heavily used in its various 
forms. These methodologies proved to be effective, although limited by their restriction 
to weak shocks and slender bodies. 

In parallel to this thrust in numerical methods, experimental investigations by 
Tijdeman 24 at the NLR in The Netherlands gave new insights into the physics of 
transonic flow, especially shock motion 2 !. Subsequent experimental investigations by 
Tijdeman et al. at NLR 42-44 and the AGARD standard aeroelastic configurations 45 ’ 46 
provided essential experimental data to be used in validation of computational fluid 
dynamics methods. In addition, unsteady aerodynamic data were also published in a 
compendium form by AGARD 47 . 

At this time, Euler and Navier-Stokes computations were still too costly, and 
the unsteady two-dimensional Euler computation performed by Magnus and 
Yoshihara 48 is noteworthy. An explicit time-marching procedure was used, with two 
Cartesian (a tine and a coarse) grids and an overlapping body-fitted grid near the 
airfoil. Another unsteady transonic computation using MacCormack’s explicit time- 
marching procedure for the Euler equations was made by Lerat and Sidfcs 4 ^. 
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A new wave of evolution in transonic flow computation came with methods for 
solution of the two-dimensional Full-Potential Equation (FPE) 50 * 52 which could handle 
arbitrary bodies, but were still limited by low transonic Mach numbers. The two- 
dimensional method of Sankar et al. 51 was subsequently extended to three- 
dimensions 4 * 5 in the USIPWING code and applied to unsteady flows past AGARD 
standard configurations 53 and a F-5 fighter wing 54 . Several improvements were also 
made in these methods, such as approximate non-reflecting far field boundary 
conditions 52 and entropy-corrected density biasing 55 * 59 . 

The two-dimensional TSD-based code LTRAN2 also benefited from 
improvements such as approximate non-reflecting far field boundary conditions 60 * 61 
and viscous corrections 62 . At the same time, research at ONERA 63 * 65 emphasized a 
strong coupling between unsteady inviscid two-dimensional TSD solutions and integral 
boundary-layer methods. Meanwhile, LTRAN2 was extended to three dimensions by 
Borland and Rizzetta 66 * 68 , in what became the XTRAN3S code, to be widely used in 
the following years 69 * 72 . Another code which evolved from this formulation was the 
CAP-TSD (Computational Aeroelasticity Program - Transonic Small Disturbance) 
code 73 * 78 , which found wide use in the research and industry communities. 


Euler methods evolved significantly with the advent of implicit schemes 14 * 79 
which allowed larger time steps. These methods were used by Levy 80 * 81 to study 
transonic buffet, and by Steger and Bailey 82 to study aileron buzz. Sankar et ai 
presented unsteady three-dimensional computations for fighter wings 83 * 84 and rotor 
blades 85 . In that method, artificial dissipation 86 is used for numerical stability. The 
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method was later extended to Navier-Stokes computations about rotor blades 6 * 87 , a 
fighter aircraft configuration 7 * 11 and a wing with an oscillatory flap 88 . An upwind 
differencing capability using Roe’s flux-difference splitting 89 * 90 has been recently 
included in that code 91 * 9 ^. Reductions in CPU times for steady transonic 1 and rotor 
blade 1 * 93 flows have been obtained by using the Navier-Stokes/Full-Potential zonal 
decomposition approach. For unsteady transonic flow, CPU times were also reduced 
with the application of the GMRES (Generalized Minimum Residual) technique 30 * 94 

A similar approach for the Euler equations was used by Guruswamy 95 . At 
NASA Ames, Holst et al. 96-98 extended Pulliam and Steger’s ARC3D code 86 to the 
solution of the thin-layer Navier-Stokes, in a new code called TNS (Transonic Navier- 
Stokes). 


Other two- and three-dimensional Euler/Navier-Stokes methods have been 
developed using upwind differencing by flux-vector splitting 99 * 100 and flux-difference 
splitting 89 * 90 . Notable examples of the former are the CFL2D and CFL3D finite 
volume codes developed at NASA Langley by Thomas et al., initially for steady 
flow 101 " 103 and later extended to unsteady flows 104 * 105 . The CFL3D code has been 
applied to and F/A-18 forebody configuration using a multiple-block approach 106-108 . 

The unstructured grid approach has received increasing attention in the past 
several years. These methods can represent virtually any complex geometry and 
adaptive meshes can be used to obtain local refinement in regions of the flow where 
gradients are larger. However, they bring additional needs for appropriate data 
structures 109 * 110 and grid refinement techniques 110 * 111 . The mesh generation itself has 


9 



been the subject of much research 110 . 112 . 1 13. At NASA Langley, Batina et al. 114 " 112 
have developed an unstructured implicit Euler finite-volume solver for unsteady 
transonic flow analysis, which has been successfully applied to a NACA 0012 airfoil 
pitching harmonically and to an ONERA M6 wing and a F/A-18 aircraft 
configuration 116 . 


To cope with the larger demands of unsteady Navier-Stokes computations for 
complex aircraft configurations, as well as to facilitate computation of moving surfaces, 
zonal structured grid approaches have been recently used 96 . 118 " 122 . Among these is the 
code ENSAERO 120 . 121 , which is in development at NASA Ames for the prediction of 
aeroelastic responses by simultaneously integrating the Euler/Navier-Stokes equations 
and the modal structural equations of motion using aeroelastically adaptive grids. These 
zonal approaches bring up the question of conservative treatment of zonal interfaces 123 . 

Problems that have received considerable attention in recent research include 
complete aircraft computations 106-108 * 116 , and delta wing oscillations 124-126 . In the 
experimental arena, a new development is NASA Langley’s Benchmark Aeroelastic 

Models Program 127 , designed to provide well documented data sets for validation of 
CFD methods. 

From the recent work on unsteady transonic flow about complex aircraft 
configurations, two trends may be identified: First, the research community is 
developing techniques for solution of the unsteady Navier-Stokes analysis over the 
complete aircraft, using appropriate approaches to deal with the complex configurations 
and large computational resources needed, mainly unstructured grids and multi-block 
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structured grids. On the other hand, there is still work to improve the transonic small- 
disturbance code CAP-TSD 128 which allows practical, fast computations that may be 
used in the design phase. This latter observation shows the gap between research and 
practical application in this area: Very complex methods being developed but a much 
simpler code being actually used. There is clearly a need to bridge the gap, by 
providing the more realistic representation of the complex Navier-Stokes methods with 
substantially reduced computational expense. The hybrid Navier-Stokes/Full-Potential 
method described in this work is an attempt at filling this gap. 



The remaining of this thesis is organized as follows: First, the mathematical and 
numerical formulation of the Navier-Stokes and Full-Potential solvers are discussed in 
Chapters II and HI, respectively. Next, the coupling between the FPE and NS solvers is 
described in Chapter IV. In Chapter V, applications of this method to a rectangular 
NACA 0012 wing in subsonic steady flow and to a F-5 wing in transonic steady and 
unsteady flow are discussed. The thesis concludes with an assessment of the hybrid 
method’s prediction capabilities and limitations and recommendations for future work. 
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CHAPTER II 


NAVIER-STOKES FORMULATION 


In the present Chapter, the mathematical formulation of the Navier-Stokes 
equations is presented and subsequently the numerical method is described 


2.1. Mathematical Fo rmulation 

The vector form of the full Reynolds-averaged, 3-D Navier-Stokes equations 
based on an arbitrary curvilinear coordinate system can be written as: 

Q r +Ef + F„ + Gf = R| + S„+Tc (2.1) 

where Q is the vector of unknown flow properties; E, F, G are the inviscid flux vectors; 
and R, S, T are the viscous flux vectors. Eq. (2.1) may be written in non-dimensional 
form, using as non-dimensionalization quantities p^ for density, for velocity, c 
(reference chord) for length, for viscosity, and p^aL for pressure, as: 

Q t +E{ + F, + G;:=-^-(R5+S,+T f ) (2.2) 
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where Re -p.^c/p.,, is the Reynolds number based on the free-stream speed of 
sound. The resulting non-dimensional flux vectors are: 


/■ > 

p 


pU 


pv 

pu 

; E = — < 

J 

puU + Z x p 

• ; F = 1 

puV + T} x p 

« pv 

pvU + Z y p 

pvV+Tt yP • 

pw 

(/ 

pwU + § t p 

J 

pwV + T] z p 

. * . 


(e+pp-lp 

y 


( e + p)V-T),p 


pW 


0 

puW + C x p 


*3 X Txx 4 *3 y Txy 4 ^ z Txz 

pvW + C y p 

r =7 

% %y'tyy + ' 

pwW + f x p 


<3x Txx 4" ^) y Tyx 4 T zz 

( e+p)W-C,p 


ZxRs+ZySs+Zjs 


0 


0 

lx Txx 4 ly Txy 4 1 z T xz 

_ 1 
; t= 7 

Cx Txx 4 Cy Txy 4 C z 't’jz 

lx Txy 4 ly T yy + Tj z fy 2 

Cx Txy + C_y 'fyy 4* Tjz 

lx Txz 4 ly Tyz 4 l z T zz 


Cx 4 Cy 4 £ z T a 

IxRs+lySs+lJs 


Cxfc + £yS5 + C z Ts 


(2.3) 


where J is the Jacobian of the transformation between Cartesian and curvilinear 
coordinates, given by: 


3 “ [ys( x $ z v-XvZs)+y v (xzzs-X(Zs)+ys(x v zz-xsz v j\ 1 


(2.4) 


U , V and W are the contravariant components of velocity, given by: 
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V = T),+T] x u+ri } v+ri,w (2.5) 

W=c,+c t u+{; r v+c.w 


The pressure p is related to the total energy e and kinetic energy by: 


p=(r- 1); 


e ~TP(« 2 + V 2 +M< 2 ) 


( 2 . 6 ) 


The shear stresses are given by: 


Txx ~f J -^i u ^ x + u v 1 lx + U;C x )^-(v^ y + Vr, T 1y + v;Cy + w^ z + Wi 1 r l z + w;Cz) 
= 4(«« ty + ' Urfy + C ,) + (Vff $ z ■ + V, I?, ■ + Vf f ,)] 

T« = Z s + Ut,T1 z + U; C m ) + {w$Cx + W v Vx+ Wff fj] 

^ == ^j( v ^ > + VT ? 77 > + Vff > )“*-(M^ ;l + M^77 x + M f f ;K + W|§ z + lV77 77 z + WfC z ) 
Tyz = A*[(v$ + V, TJ g + Vf f 2 ) + (w$ £ y + Wtj *7, + W? £,)] 

+ w u + wcf,)' - 1(«$ {,■ + n x + U; c x + V(C,- + Vt, r]y + V;Cy) 


J (2.7) 


Ta 


and 


Rs = ui:„+vr v +wT a +-^-^^ x d i a 2 +V,d 1 a 2 +C x d ( ‘‘ z ) 

S5 = UT v + vt„ + wt ) , + — (g, dt a 2 + % «?„ a 2 + ?, «? f a 2 ) (2.8) 

r5 = «r»+VT,i+tVT a +^— -jy (£.<?{ a 2 + r l,d 1 a 2 +^,d(a i ) 
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where Pr = fic P /k is the Prandtl number and a is the speed of sound. The notation 
d$ a 2 is a short form for d{a 2 )/d^ t 

In turbulent flows, the molecular viscosity fi appearing in Eqs. (2.7) and (2.8) is 
replaced by \i +/z r , and the quantity /z/Pr in Eq. (2.8) is replaced by /z/Pr +H T /PT Ti 
where fi T is an eddy viscosity and Prr is the turbulent Prandtl number. 


2.2. Numerical Formulation 

In this section the finite-difference numerical formulation of the Navier-Stokes 
equations (2.2) is discussed. First, the finite-difference discretization of the derivatives 
is described. Next, the linearization of the resulting non-linear system of equations and 
its approximate factorization into two block-tridiagonal systems of equations are 
discussed. Finally the numerical implementation is described. 


2.2.1. Discretization 


The time derivative, Q x , of equation (2.2) is approximated using two-point 
backward difference at the new time level 'n+T: 




(2.9) 
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where n refers to the time level at which all quantities are known, and ‘n+1’ is the 
new time level. All spatial derivatives are approximated by standard second-order 
central differences and are represented by the difference operators 5 , e.g.: 

(Ef)„, = + 0( Af ) = + 0 ( A?) (2.10) 

Note that the choice of = A 77 = A£ = 1 in the computational space is made 
for convenience. 


The streamwise and normal derivatives, and G^, are evaluated implicitly at 
the new time level 'n+1'. The spanwise derivative, F^, is evaluated explicidy at the old 
time level 'n', but uses the ‘n+1’ values as soon as they become available. Thus, the left 
side of the discretized form of Eq. (2.2) becomes 


| BgM-EEfo , KW-gTfejt G"!w-G'u-i 
At 2 2 2 


( 2 . 11 ) 


This semi-explicit treatment of the spanwise derivative enables the scheme to 
solve implicidy for AQ n+1 at all points at a given spanwise station at a time. To 
eliminate any dependency the solution may have on the sweeping direction, the solver 
reverses the direction of spanwise sweeping with every sweep, i.e. for every other 
sweep, Eq. (2. 11) is replaced by: 


Qw~Q*m 1 EfflM-E; ilu | F&u-FL-u Gfiki-GS! 
At 2 2 2 


( 2 . 12 ) 
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The viscous terms R^, and are evaluated explicitly, using half-point 

central differences denoted here by the difference operator 8 , so that the computational 

stencil for the stress terms uses only three nodes in each of the three directions. For 
example, in the computation of R^, the term IfAu ^ appears and is discretized as 

follows: 


= 8(?J*t)+0{. Af )= W** +0 (At 2 ) 


(F 2 iA Mi+l./.fc Mi./,* ( e2..\ Ui.j.k ~ Ui-lj.k 

^ Ly 2 ,j.k A g 

i^‘^\j,k + »^)i+l,/,t (ui+l,J,k - Mi,/,*) (^)i,;,* + (^)i-u.k (uuj, - Ui-lJ,k) 

2 M 2 






(2.13) 


2 + ‘4 + , /,*](««■>.* - «U.t) “ + (l ^),_i “ ««,/.*)} 


Explicit treatment of the stress terms still permits the use of large time steps 
since the Reynolds numbers of interest here are fairly large. 

With the above described time and space discretizations, the discretized form of 
Eq. (2.2) becomes: 


AQ” +1 + At(5 4 E" +1 + 5„F"'" +1 + <S ? G” +1 ) = |%R-" +1 + s„ S”-" +I + SjT"-”* 1 ) (2.14) 


Note that all viscous terms include T| -derivatives, for which known values at the 
new time level 'n+T are used, hence the notation ^R n, " +1 . 


17 



2.2.2. Linearization 


The time and space discretizations described above lead to a system of non- 
linear, block penta-diagonal matrix equations for the unknown AQ" +1 = Q n+1 - Q n , Eq. 

(2.14), since the convection fluxes E, F, G are non-linear functions of the vector of 
unknown flow properties Q. Equation (2.14) is then linearized using the Jacobean 
matrices A = dE/dQ and C = 9G/dQ, given by: 


( W 



kt 

k x <t> 2 -u0 

ky<p 2 -vd 

k z (j> 2 -wd 
1 9(<p 2 -b) 


k x 

ky 

kz 

0 

®-kx7 2 u 

kyM kx Y\ v 

kzU-kxYiW 

kxY\ 

kxV-kyYiU 

®-kyY 2 v 

kzV-kyYiW 

kyY\ 

kxW kzY\U 

kyW-kzYiV 

®~kz Y 2 w 

kzYi 

kxb-YiuO 

kyb~YivO 

kzb-Y^wQ 

kt+Y& 


(2.15) 


where: 



for A 
for C 


® = k t + 0 


^ 2 = -^— ( m 2 + v 2 + w 2 ) ; 9 = k x u+kyV+kzW 

Yi=r-i ; r 2 =y -2 ; (2.16) 

p 


The linearization is obtained as follows: 



(2.17) 
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and similarly for G n+1 . Applying these linearized flux vectors to Eq. (2.14) yields: 
[l + A" + 5 {C")]aQ” +1 = RHS”'“ +1 

= -At(5 { E“+ <5,F"'" +1 + S(G") + ~(S( R“" +1 + s v S"^‘ + S ( T"’ +1 ) * 

1AV ' 


where I is the identity matrix. 


2.2.3. Approximate Factorization 

Eq. (2.18) is a system of linear, block penta-diagonal matrix equations, which is 
considerably expensive to solve. The approach used here is to employ the approximate 
factorization of Beam and Warming I 4 : 

[1 + At(<^ A* + 8; C n )] AQ n+1 = [I + A T$s A n ][l + At 8; C n ] AQ B+1 + 0(At 3 ) (2.19) 

which allows the system of equations (2.18) to be written as: 

[I + A t8s A B ][l + At 8; C n ] AQ B+1 = RHS n ’ B+1 (2.20) 

Note that there is no loss of temporal accuracy, because the error incurred due to 
the approximate factorization is of order 0 (At 3 ). The system of equations (2.20) may 

now be solved in two steps, each involving only a block- tridiagonal system of 
equations: 
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[I + At£| A n ] AQ*" +1 = RHS n,n+1 
[I + At^C*] AQ B+1 = AQ* n+1 


( 2 * 21 ) 


2.2.4. Diagonal Form 

The computational work required to solve the systems of equations (2.21) may 
be reduced by employing the diagonal algorithm of Pulliam and Chaussee 13 . Since the 
flux Jacobian matrices A and C have a complete set of eigenvalues and a corresponding 
set of distinct eigenvectors, the similarity transformations may be used to diagonalize 
A andC: 


A = T*A*T f ; C = T f AfTf 1 


( 2 . 22 ) 


where the diagonal matrices A^ and A^ may be concisely expressed as: 


A$ = diagjtf, £/,£/,£/ + a^,U - a->jX\ 

A f = dia g[\V, W, W, W + a^A*, W - a^ 6 ] 


(2.23) 


where A\-t; x + t; 2 + t; 2 t and Ae - f 2 X + f 2 + C 2 . The eigenvector matrices are given by: 
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T* = 

kx 

L» 

kx v + kiP 
kx w -k y P 

[kx<t> 2 /7i + 


ky 

kyU-k,P 

k, v 

kyW + LP 
k y <t> 2 /ri+ 


k z 

kz u + kyp 

kzV-kxP 

kz w 

[k,ftYi + 


p(i.v -*,")] p{k x w-k,u)] p(kyU-k x v ) ] 


a 

a(u+k x a) 

a { v+ k, a ) 

a(w+k t a) 


a 


L x. 


4 - dO 


a 


a 

a(u-k x a) 

«(v-k,a) 

a ( w ~ic, a ) 

f+a 2 

L r, 


-aO 


(2.24) 


and 


T Ji l = 


1 

X 

‘"s 

1 

M 

i 

kx u S 

kyf-P' l {L w ~k, u ) 

kyUg-k.P' 1 

k,f-P~ l [kyU-k x v ) 

k*“g + kyP~' 

p{f-ad) 

P{k x a-Y i«) 

p(<p 2 + a9) 

-Pik.a+Yiu) 


k^g + hP* 

kxWg~kyP~' 

~k x 8 

kyVg 

kyWg+k x P~ l 

-fCyS 

k.vg-LP' 1 

k.wg 

~kz8 

P(ky a ~Yiv) 

P{k, a ~r, w ) 

PYi 

~P(ic y a +Yl v) 

~P{k, a + Yi w ) 

PYi 


(2.25) 


where 


0 


a -J2a ’ ^ V2 pa 

k.u+kyV+k.w ; f = (l-<l> 2 /a 2 ) 


'X,y,z 


V kl + k 2 + k z 


; g=(r-i)/a 2 


(2.26) 


Applying the similarity transformations (2.22) to Eq. (2.20) and using the 
identities l = T k T k 1 yields: 
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[(T4T{ 1 )'+AT5{(T{A ? TJ 1 )"][(TfTf‘)’ + Ar ( 5 f (T f A ! :Tj 1 )"]AQ" +1 = RHS“'” +l (2.27) 


The fundamental simplification of the diagonal algorithm is obtained by moving 
T$ and T$ outside the difference operators 8$ and 8;, respectively: 

Tj[l+ A t 5 ? A|]N"[l + At<5 ? A^T^’AQ**' = RHS"-" +1 (2.28) 

where N = T| l T{. This simplification introduces and error because T$ and T; are 
functions of and cannot be arbitrarily brought in or out the derivatives. It is 

however believed 13 that the errors introduced are of order 0( At 2 ) and have the effect 

of making the scheme first order accurate in time, the same order as the previous 
approximations. 

The solution of Eq. (2.28) still involves two block-tridiagonal systems, but now 
the blocks are diagonal matrices. The solution of Eq. (2.28) is obtained through the 
following steps: 


[l+ Arg s A{]AQ** +1 = (T{ 1 )”rHS" - " +1 

[i + A t 8; A{] AQ**" +1 = ( N -‘)"AQ*“ m (2.29) 

AQ" +1 = T£AQ"“ +1 

Note that the interim variables AQ and AQ may be stored in the same memory 
locations as AQ , conserving the available memory. 
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2.2.5. Numerical Dissipation 


The use of standard central differences to approximate the spatial derivatives 
can give rise to the growth of high frequency errors in the numerical solution with time. 
To control this growth, a set of 2nd/4th order non-linear, spectral radius based, explicit 
artificial dissipation terms are added to the discretized equations. A second order 
implicit dissipation is used to help the overall numerical stability of the scheme 86 * 129 . 
The systems of equations in (2.29) are modified as follows: 

[l + - At£/V$ 0i At /]AQ*' l+1 = (Ti 1 )"(RHS" ,n+1 - D?*" +1 ) 

(2 30) 

[1 + At^ AJ - Are/ Vf 03 Af/] AQ** rt+1 = (n-^AQ^ 1 

where e 7 is the coefficient of implicit numerical dissipation, V$ and A$ denote 
backward and forward difference operators, respectively, <p l and 0 3 are defined as: 


2 ' 1 ' 


14- 


max 


V+d'sfA 4 W+a^[A6 
U+ajAS U+ajAt ) i+l 


1+ max 


U+a^Ax V A4 
W+d'sfXe W+a^fAs 


k+— 
2 


(2.31) 


and D" ,n+1 is the explicit fourth order dissipation, given by: 


D"’" +1 = At e £ [v 4 & (a 5 V { A 4 ) + 0 2 ( V,A,,) 2 + V 4 0 3 (A { V 4 Aj)]/ Q“ (2.32) 


where ee is the coefficient of explicit numerical dissipation and 
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1 + 


f 


max 


v 


U + a'\[Ai W + a^fA6 
V + a^'V + aJT 4 


(2.33) 


For points adjacent to the computational boundaries, second-order explicit 
dissipation is used instead of Eq. (2.32). 

The above described fourth order dissipation may lead to “wiggles” near shocks. 
To avoid this problem, a switching function based on the second normalized 
streamwise derivative of pressure 


Pm-ZPi + Pm 

\Pi + i~2Pi + Pi-i 


is used to replace the fourth order dissipation with second order dissipation near 
shocks 130 * 131 . 


2.2.6. Turbulence Mode! 

A slightly modified version of the Baldwin-Lomax (B-L) algebraic turbulence 
model 132 is used, where the maximum shear stress is used instead of the wall shear 
stress because in the vicinity of separation points, the shear stress values approach zero 
at the wall. 

In this model, two layers are considered; in the inner layer, \i T is given by: 
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(^r)inn«r layer 


( 2 . 34 ) 


where |fi)| is the mean vorticity, given by: 


(dw_dv) ^ du 
dy dz ) 


dv du f 
dz dx ) ^ dx dy j 


(2.35) 


and £ m is the mixing length, given by: 


£ m = Kd[ 


(2.36) 


where *c=0.41 is the von Karman constant, d is the distance from the wall, 
A+ = 26.0 is the van Driest constant, and 


d+ = d 



(2.37) 


The modification with respect to the original Balwin-Lomax model is apparent in 
Eq. (2.37), where Tma* is used instead of Twaii. 

In the outer layer, fi T is given by: 

(^r) outcr layc. = KcpCiFwFk (2.38) 

where Kc = 0.0168 is Clauser's constant, c\ - 1.6 is an empirical constant, F w is given 
by: 
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F w = min dm^F maXi 0.25 


(239) 



where 


F 


max 



(2.40) 


Ud if = ma x(Vm 2 + v 2 + w 2 ) - min(Vw 2 +v 2 + w 2 ) 


(2.41) 


and (i m is the distance from the wall where F max occurs. Also in Eq. (2.38) Ft is 
given by: 


F k = 


1 


1 + 5.5 



(2.42) 


The switch between inner and outer zones occurs at the distance d c , defined as the 
smallest distance from the wall for which (^r) iBralayCT = (/ir) ou „ tUyct . i-e., the values from 

Eqs. (2.34) and (2.38) are the same. 


2.2.7. Numerical Boundary Conditions 


The formulation described above must be complemented by appropriate 
boundary conditions to be specified along the solid surface, Full-Potential/Navier- 
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Stokes interface, at the wing root, and far field boundaries located outboard of the wing 
tip and downstream or outflow boundary beyond the wing trailing edge. The boundary 
conditions to be applied at the Full-Potential/Navier-Stokes interface are discussed in 
Chapter IV. The remaining boundary conditions are discussed next. 

2.2.7 .1. Solid Surface Boundary 


The solid surface corresponds to the plane k = 1. The unknown vector in Eq. 
(2.29), AQ includes values from k = 2 to k — KMATCH . At the end of each iteration, 
the new values of Q, w>1 are computed as follows: Density and pressure are computed 

from the assumption that their normal derivative at the solid surface is zero, 
dpj dn = dp/ <9/i = 0. This is approximately satisfied on near-orthogonal grids as: 


O — ^ Pi,;, 2 Pi,/, 3 


Pw 3 


(2.43) 


The velocities at the surface are computed from the no-slip condition, i.e.: 

Mi,/\i = (xx) iJtl Vij,i = (y r ). A1 wij, i = {zx) iJtl (2.44) 


2.2.7.2. Wing Roof Boundary 

The wing root corresponds to the plane j = 1. The unknown vector in Eq. 
(2.29), AQ includes values from j = 2 to j = JMAX-l. The values of Q- u , i.e., at 
the root, are not updated; when computing the residual RHS' , *" +1 at the j = 2 cell, the 
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fluxes at j-1 are computed with the symmetry condition that the contravariant 
velocity normal to the boundary vanishes, i.e., V = 0. The pressure values at j = 1 are 
computed using zeroth order extrapolation, so that p n k = p. 2 k . 

2.2.7.3. Far Field Boundaries 

The downstream (* = 1 and i = IMAX) and outboard {j-JMAX) boundaries 
are treated in the same way. The velocity normal to the boundary is computed. Then, 
the boundary conditions are imposed depending on whether it is an inflow or outflow 
and whether it is subsonic or supersonic: 

• Supersonic outflow: all variables are extrapolated from the interior of the 
domain; 

• Subsonic outflow: the pressure is fixed to be the free-stream value and the 
other variables are extrapolated; 

• Subsonic inflow: the density is extrapolated from the interior of the domain 
and the other variables are fixed from the free-stream; 

• Supersonic inflow: all variables are fixed to be the free-stream values. 
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CHAPTER III 


FULL-POTENTIAL FORMULATION 



In the present Chapter, the mathematical formulation of the Full-Potential 
equation is presented and subsequently the numerical method is described The Full- 
Potential solver used in the present work was developed by Sankar et al. 4 ^,133 


3.1. Mathematical Formulation 


The 3-D unsteady compressible potential flow equation, in a body-fitted 
coordinate system, may be written in a strong conservation form as: 



(3.1) 


where p is density and U , V and W are the contravariant components of velocity, 
given by: 
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v= T h+n x u+iiyV+v,w 
W = t,+Z x u+Z s v+Z I w 


( 3 . 2 ) 


End J is the Jacobian of the tran s f orm a tion between Cartesian End curvilinesr 
coordinates, given by: 

J = [y(( x ^’i-x^;)+y^x i z ( -x ( Zi)+y ! (x n z f -X 42 n)}~' ( 3 . 3 ) 

In the present formulation, the full potential is denoted by tp and the 
perturbation potential is denoted by (p , i.e.: 

v = </> y =v M + (p y (3.4) 

w - (j) z - + (p z 

It should be noted that the contravariant components of velocity can be 
expressed in terms of the derivatives of <j> by substituting (3.4) into (3.2), which yields: 

U- £,+ Al0£ + A2 0,, + A3 0f 

v = T) t + Ai <Ps + Aa (p v + As ( 3 . 5 ) 

^ = + + A5<P v + A6<I>s 


where 
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A2=Z t V, + S y ri y + S I ll, 

A^it+rtl+rj 1 , 

As = ih£ x +il y C y +Ti,£, 
Ae = C] + £ + C 


In addition to the differential equation (3.1), an additional relation is needed to 
express the density in terms of the velocity potential and its derivatives (i.e., the 
velocity components). This additional relation is the isentropic gas law 

1 

7-1 

(3.7) 



where a is the speed of sound, given by the energy equation: 


a 1 , a , m 2 + v 2 + w 2 
. + <*>,+ 


7-1 


at V 2 
— — + — 
y- 1 2 


(3.8) 


Note that the derivative (p t may be expressed in terms of the derivatives with 
respect to the transformed variables as: 

9t = 9t + Vs Zt + Vr, V t +V; C t (3.9) 


Using Eqs. (3.7) and (3.8), Eq. (3.1) may be written 4 as a second order 
hyperbolic partial differential equation for <p : 
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( 3 . 10 ) 


^j[<Prr+ + V< Pv + "%] = + Q 

where Q is a source term associated with the rate of grid deformation 133 , which 
vanishes for rigid grids and may be neglected for mildly deformed grids, which is the 
case considered here. 


3.2. Numerical Formulation 

In this section the finite-difference numerical formulation of the Full-Potential 
Equation (3.10) is discussed. The spatial flux-like terms appearing on the right hand 
side of Eq. (3.10) are discretized using standard central differences, which result in 
formal second-order accuracy in space. The mixed time-space terms appearing on the 
left hand side are discretized using two-point upwind differences. The temporal 
derivatives are discretized using two-point backward differences. These discretizations 
are described in more detail below. For convenience, the mesh spacing A§ , Atj and 
Af are set equal to unity in the computational domain. 


3.2.1. Discretization 

At a given time level w, the disturbance velocity potential (p and its temporal 
derivative (p x are known, and consequently all velocity components, speed of sound 

and density are also known. Eq. (3.10) is a partial differential equation for <p with 
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nonlinear coefficients. To circumvent the nonlinearities, the coefficients p, J , U , 
V and W appearing on the left side, and the density p appearing on the right side of 
equation (3.10) are computed at the time level n. The remaining quantities in (3.10) are 
kept at the new time level n + 1. In the process of evaluating the contravariant 
velocities U , V and W , two-point central differences are used to evaluate the 
derivatives of (p and the transformation metrics at the grid points and locations mid 
distance between the grid points. 


The temporal derivatives on the left hand side of Eq. (3.10) are discretized using 
two-point backward finite-difference operators, while the mixed time-space terms 
appearing on the left hand side are discretized using two-point upwind differences. In 
this respect, the left hand side of Eq. (3.10) is expressed as follows: 


+ U"S ( 8 r <p + V'S^ip + ( 3 . 11 ) 


For example, at a typical grid node (i,j,k), the first term inside the square 
brackets of Eq. (3.11) is expressed as 


S g - 9** ~ 2( P n ± <P n ~ l _ &<p n+l ~ Aff" 
(At) 2 (At) 2 


(3.12) 


In the previous expression, A (p represents the change in the solution in two 
consecutive time steps, i.e., A <p nn = <p n+1 - <p n , and A <p n = q> n - <p n ~\ The mixed space- 
time derivatives appearing in Eq. (3.11) are discretized using upwind-differencing for 
the spatial derivative, and two-point backward-differencing for the temporal derivative. 
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For example, in evaluating the second, third and fourth terms, the following 
expressions are used: 


wi,-£± M 


+ u-\u\ 

rAC‘,- A cv 


s 2At 

W8.8 a - W+ \W 

A§ J 

irACi-A^., 

2Ar 
~ + w-\\ 

L 

Yl 

[a^L.-ac 

i ■ 
k_ 

2At 

L ^ 

AC*-AC-ul 

2Ax 

. v-M 

‘a 

Af 

cc-ac.*1 

_ 

" **“ 2At 

Ar\ 

T ' 

2Ar 

At] j 



(3.13) 


The flux-like terms appearing on the right hand side of Eq. (3.10) are evaluated 
using two-point central-difference formulas, i.e.. 



The density p in Eq. (3.14) is computed at the time level n, while the 
contravariant components of velocity are computed using mixed information from time 
level n and the new time level n + 1, in order to reduce the number of diagonals in the 
final matrix of coefficients. Recalling Eq. (3.5), the contravariant components are 
evaluated as follows: 
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U’« = l + Mf' + A2<l> n v + A 3 <l>} 

= Z, + Ai(<t>" f + A<p"*') +A 2 <t>v + A 3 ipz = U'+A l A 1 

v = n.+Ai^+A^'+Asr; 

= V, + Arf] + A t (r, + Acp” 1 ) + A 5 r ( = V" + A 4 A(p'« 

= C+Aa^ + ^ + A^* 1 
= C, + A 3 <t>l + As K + A«(0f + Ap” +I ) = iy" + A 6 A^j +l 


3.2.2. Density Biasing 


In order to maintain numerical stability in regions of supersonic flow, the 
numerical formulation must be constructed in such a way that it is consistent with the 
physical domain of dependence. For that purpose, the artificial compressibility 
method 134 is used. Here, the density values p that appear in (pU/J) on the right side 

of equation (3. 10) are biased in the direction of the flow using a procedure suggested by 
Hafez, Whitlow and Osher 57 . First, a function F is defined as: 

F = pq if Af > 1 

C * -c (3.16) 

F=p q if M^l 

where the superscript * refers to sonic conditions. 


Then the biased density is defined as: 


Pi „ = P i 
l+ 2' 1,k 


*V’* 


-F 


l ~2' hk 


Qij.k 


(3.17) 
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where q is the flow speed and M is the Mach number. It is clear from Eqs. (3.16,17) 
that the biased density reduces to the local density in subsonic flow regions. 


3.2.3. Strongly Implicit Procedure 


When the above discretizations are employed, at each grid point a linear 
equation results for the quantity A(p n+1 = (p n+l - <p n y namely: 


+e?j, k A(p*;lj tk +f*j ik Atf* l +Uk + g-j tk A(p-j Ml = Rfjjt 


(3.18) 


where the coefficients a?j, ki b?j, k , Cij, kl d?j, k , e*j, k , f* jik , g* j k and Rfj tk are functions 
of the transformation metrics, the contravariant velocities, the density p, the speed of 
sound a, and the time step At. Application of Eq. (3.18) at the grid points result in a 
sparse pentadiagonal matrix system which may be expressed as: 


[M]{Acpr l ={RY 


(3.19) 


A lower-upper (LU) approximate factorization scheme, originally devised by 
Stone 135 , and applied to transonic flows by Sankar, Malone and Tassa 4 , is employed to 
solve the system of equations (3.19) efficiently. In Stone's strongly implicit procedure 
(SIP), the matrix [M] is approximately factored as the product of two sparse lower 
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([L]) and upper ([£/]) matrices each having four diagonals. Eq. (3.19) can thus be 
expressed as: 

[L][U]{A(pY +1 ={R} H (3.20) 

where the elements of matrices [L] and [U] are recursively related to the coefficients 
of the matrix [M]4>133. The solution to Eq. (3.20) is then obtained using a two-step 
procedure where we first solve for a temporary solution vector i.e.: 

[L]{A(p'} = {R} n (3.21) 


and next solve for {A<p}" +1 : 


[C/]{A«,r={A^} 


(3.22) 


It should be noted that the above approximate factorization procedure is 
applicable to both quasi-steady as well as unsteady flow field solutions. In the former 
case, the temporal derivatives of the potential function are set equal to zero and the SIP 
can be regarded as an iterative relaxation procedure. In the latter case, the SIP is 
regarded as a one-step non-iterative time-accurate marching procedure. 
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3.2.4. Numerical Boundary Conditions 


The formulation described above must be complemented by appropriate 
boundary conditions to be specified along the FuH-Potential/Navier-Stokes interface, at 
the wing root, and far field boundaries located outboard of the wing tip, outer 
boundary, downstream or outflow boundary beyond the wing trailing edge. The 
boundary conditions to be applied at the FuU-Potential/Navier-Stokes interface are 
discussed in Chapter IV. The wing root and far field boundary conditions are discussed 
next. 


3.2.4.I. Wing Root Boundary 

The governing equations (3.10) are applied on the j = 2 cell as at other interior 
points. The computational plane j = 1 corresponds to the wing root. At this plane of 
symmetry the contravariant component of velocity V should be zero. As an 
approximation, this condition is enforced at j-2. Using the expression for V in Eq. 
(3.5), and the condition V = 0, the derivative <p v may be obtained as 

tlMl (3.23) 

Aa 

After 0, j is found, the derivative of the perturbation velocity potential, <p v is 
easily computed and the perturbation velocity potential at j = 1 is computed from 

%.u = - 2(Ar,)[*,„]. 2i (3.24) 
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3.2.4.2. Far Field Boundary 


Since the flow is considered to be uniform at a large distance from the wing, the 
disturbance velocity potential function <p is usually set to zero on all far field 
boundaries. This condition also implies that the flow velocities in the planes containing 
these boundaries assume free stream values. 

In practice, the computational domain is bounded and the assumption of zero 
disturbance potential at distances not very far from solid surfaces causes the acoustic 
waves that carry the perturbation information to be reflected at the outer computational 
boundary. These reflected waves contaminate the solution and delay convergence 136 . 

In order to minimize this adverse effect, the far-field boundary conditions 
derived by Shankar et al.^2 are used at the outer boundary. In this approach, the 
Riemann invariant R that corresponds to positive characteristics with respect to the 
inward normal to the boundary is specified. The Riemann invariant may be expressed 
as: 


R = - 


W 

V^6 



(3.25) 


The actual implementation of this approximate non-reflecting boundary 
condition is earned out as follows: Let R™ be the Riemann invariant corresponding to 
the undisturbed flow field, i.e.: 


— 


W. , 2 

-f. doo 

4 m r-i 


(3.26) 
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The internal solution in the computational domain at a given time level n gives 
a Riemann invariant R&kmax- i at the point (i,j,KMAX - 1). Variations in the values of 
(p and its derivatives will introduce a variation SR with respect to R?j,kmax- j . The 
objective here is to set the values of the perturbation potential <p at the outer boundary 
so that at each time step the variation that corresponds to these new values of <p is 
equal to the difference between /?« and RI},kmax~\ , i.e.: 

8Ri,j,KMAX- 1 — ~ Ri,j,KMAX~\ (3.27) 


so that the Riemann invariant in the computational domain approaches the far field 
value. By employing a variational calculation on the expressions of W , Eq. (3.5), and 
a , Eq. (3.8), and neglecting variations in the tangential derivatives of the potential, 

A A A 

5(p^ and 8(p^ the following expression may be obtained for SR: 


SR - -4M Sq)^ ~ 5(p^ + 5(p)j 


(3.28) 


The following difference approximations are used: 


(*^L .KUM~ A<P ‘^ KMMC 


„ ^ViJ.KMAX ^ ( PiJ,KMAX 

t/i,j,KMAX fa? 


(3.29) 


It should be noted that the change in notation in Eq. (3.29) from S(p to A<p 
implicitly assumes that the variations in the derivatives are due to the changes in (p as 
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the solution marches in time, since A (p n+l = (p n+l -(p n . By applying Eqs. (3.28) and 
(3.29) in Eq. (3.27), the following expression may be obtained for A 


^ViJ.KMAX 


1 + A + ajAl) + Ar i W ~ + a ^) ~ aAt i R - ~ 


(3.30) 


Eq. (3.30) is used at each time step to update the values of <p at the outer 
boundary. Similar approaches may be used for the downstream (/ = 1, i = /MAX) and 
outboard (j = JMAX ) boundaries. However, in the present work the application of the 
above described boundary conditions only on the k = KMAX boundary was sufficient 
to give good convergence characteristics. 



CHAPTER IV 


NAVIER-STOKES/FULL POTENTIAL COUPLING 



A typical partitioning of the domain into an inner zone and an outer zone is 
illustrated in Fig. 4.1. The plane k — KMATCH corresponds to the interface between 
the inner zone and the outer zone. The Navier-Stokes solver is applied at all planes up 
to k = KMATCH. The FPE solver is applied between the planes k = KMATCH and the 
outer boundary k-plane. Therefore, the two zones actually overlap, which allows 
specification of boundary conditions at the interface without extrapolation. 
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Fig. 4.1: Partitioning of Computational Domain into Inner and Outer Zones. 


4.2. Hrigjpal Yiscous-Inviscid Interface Boun dary Conditions 

In this section, the interface boundary conditions used in the original hybrid 
scheme 1 are described. The Navier-Stokes solver is applied up to the location 
k = KMATCH . This solver requires the flow properties (density, velocity, pressure) at 
the plane KMATCH + 1. These values are obtained through the numerical 
differentiation of the velocity potential, and the application of the isentropic energy 
equation: 
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(4.1) 





(4.2) 


The full potential solver is applied down to and including the plane 
k- KMATCH . This solver requires the velocity potential at the plane 
k = KMATCH - 1 . These values are obtained by matching the normal component of 

velocity v n at the plane k = KMATCH from the potential flow and the viscous flow. In 
terms of the perturbation potential, this is equivalent to matching at k = KMATCH 


and using the difference formula: 



i,j,KMATCH 


$ i ,j, KMATCH +1 ij, KMATCH -l 

2 


(4.3) 


In order to obtain (p^ from the Navier-Stokes solution, the expressions for the 

contravariant component of velocity W in terms of the primitive variables and 
perturbation potential are used: 

W^C'H'U+CyV+C.w (4.4) 

w = C « + C,w- + ? , V- + C,w. + A 3 <Pf + As <p v + As <P; (4.5) 
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Eq. (4,4) is used to evaluate W in terms of the primitive variables, obtained 
from the Navier-Stokes solution at k = KMATCH. Using the value of W from Eq. 
(4.4) and (p^ and (p v from the numerical differentiation of the perturbation potential 
(p at the previous time level V, <p ; is obtained from Eq. (4.5) and finally used in Eq. 
(4.3) to obtain <p at k = KMATCH - 1 . 


4.3. Improved Viscous-Inviscid Interfac e Boundary Conditions 

Previous applications of the hybrid NS/FPE solver to an iced wing 
configuration 1^7 with the above interface boundary conditions showed an oscillatory 
behavior in convergence histories that indicated false reflections from the Navier- 
Stokes/Full-Potential interface when the boundary conditions where implemented as 
above. Similar numerical phenomena were observed in the past with respect to far-field 
boundary conditions: Acoustic waves traveled from the solid surface to the outer 
boundary and were reflected back to contaminate the solution and delay 
convergence! 36 . The spurious waves responsible for the oscillatory convergence 
behavior need to be eliminated. In unsteady flows this is even more important since 
these spurious waves will compromise the time accuracy of the solution. 

In the past several years, research has been underway 52 » 60 » 6 l» 1 38-150 10 develop 
and apply non-reflecting far field boundary conditions, which accelerate convergence 
to steady-state and in some cases improve time accuracy. These boundary conditions 
would not be directly applicable to the viscous/inviscid interface discussed here 
because perturbations in one zone must be transmitted to the other zone, as illustrated 
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in Fig. 4.2. The disturbances in the inner region should contribute to the outgoing 
waves only, while the disturbances in the outer region should contribute to the ingoing 
waves only, so that there is no reflection at the interface. 



Full-Potential node 


Outgoing waves 



Incoming waves 


k = KMATCH 
interface 



Navier-Stokes node 


Fig. 4.2: Waves Contributing to Fluxes at k - KMATCH Interface. 


4.3.1. Interface Boundary Conditions for the Navier-Stokes Solver 

Following a development analogous to Giles’ derivation of approximate non- 
reflecting boundary conditions 148 , a set of characteristics normal to a f ^constant 
surface was obtained as follows: First, the vector form of the 3-D Euler equations based 
on an arbitrary curvilinear coordinate system can be written as: 

Q r +E$ + F„ + G f = 0 (4.6) 
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where Q is the vector of unknown flow properties; E, F, G are the inviscid flux vectors. 
These vectors are given in Eq. (2.3) and are repeated here for convenience; 


p 


pU 

pu 

; E = ~ 

/ 

puU + Z x p 

pv 

pvu+4 y p . 

pw 

J 

pwu+4,p 

. * . 


(e+p)U-4 t p 


pv 


pW 

puV +TJ X P 

; G = — ’ 

J 

puw+4 xP 

pvV + r\ y p 

pvW + £ y p ■ 

pwV + Ij'P 

J 

pww+4 s p 

(e + p)F-?7,p 


(e + p)W-t,p 


(4.7) 


Next, small perturbations on the primitive variables are considered. Let the 
vector of perturbations on the primitive variables be denoted by <5q, where: 

5q = [8p Su 8v 5w 8pf (4.8) 


Then the small perturbation form of the Euler equations (4.6) is given by: 



where A, B , and C are Jacobian matrices for the small perturbations. Here, the 
purpose is to construct a set of characteristics normal to a £ =constant surface, therefore 
tangential variations are neglected, and consequently Eq. (4.9) is reduced to: 

<5q T +C<5q ? = 0 (4.10) 

where 

W pC x Pty PC Z 0 

o w 0 0 

p 

r 

o o w o 

p 

o o o w ^ 

p 

o 7PC X wCy wC z w 

The matrix C has five eigenvalues: Xi = X 2 = Xi = W, X* = W -a^As, and 
Xs = W + a-\[jU , where As is given in Eq. (3.6). The five characteristics corresponding 
to the hyperbolic system (4.10) are constructed by applying the similarity 
transformation: 

C = t { A?ff 1 (4.12) 

Applying (4.12) to (4.10) and left-multiplying the result with f^ 1 the 
characteristic equations are obtained: 
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Cr + AfCf = 0 


(4.13) 


where c = fj'5q , i.e 


' 

Cl 


0 tf. 

J,C. 

Cx + Cy 

0 

sp 


Ae 

A6 

As 



C2 


1 

o 

Cl+C 2 , 

Zytx 

0 

Su 



A6 

A6 

Ae 


c 3 

► = 

1 0 

0 

0 

1 

a 2 

Sv 

c 4 


Q P a Cx 

2 Ja<, 

pa£ y 

2VaI 

pa£. 

J_ 

2 

Sw 

r% 

Cs 


n 

paCy 

P a C, 

1 

Sp 

s 




2 4a 6 

2 

* - 


(4.14) 


The integration of Eq. (4.13) to obtain the characteristics c\ to cs is performed 
according to the signs of the corresponding eigenvalues, for example if yt 4 > 0 then c 4 
is computed using information from the inner zone, otherwise it is computed using 

information from the outer zone. This corresponds to the eigenvalue splitting 
Af = AJ + Af and corresponding characteristic splitting c = c + + c". Then Eq. (4.13) 


becomes: 


ct+A|c£ = 0 (4.15a) 

c^+AfCj = 0 (4.15b) 

The integration of Eq. (4.15a) is discussed next For convenience, one scalar 
equation is treated. Let ct denote the characteristic at the k plane being calculated (the 
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subscripts ij are dropped to shorten the notation). The discretization of the scalar form 
of Eq. (4.15a) is made as: 


/ , \new / yu / + \new / . \NS 

At Af 


(4.16) 


where {ct-i) NS is the characteristic given by the Navier-Stokes solver at the k - 1 plane. 

Considering that the linearization (4.9) is performed about the previous time step, then 
(rf)* = 0 and Eq. (4.16) can be solved for (c*) new as: 


(cty 


A* At 
A + At+1 



(4.17) 


The integration of Eq. (4.15b) is performed in an analogous fashion. The 
discretization of the scalar form of Eq. (4.15b) is made as: 


/ _\new / _\old / _ \FP / _\ncw 

""(g*) f -(c*) q 

At Af 


(4.18) 


where (c*+i) FP is the characteristic given by the Full-Potential solver at the k + 1 plane. 
Again, considering (c*)° ld = 0, Eq. (4.18) can be solved for (c~ k f™ as: 

For steady flows, it has been found that slightly higher convergence rates may 
be obtained by taking the limit as At -> °o in Eqs. (4. 17) and (4. 19), i.e.: 
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(420a) 

(4.20b) 


(ctr = (cU) NS 

(c*r=(c- M f 

For unsteady flows, numerical experimentation has shown that using the 
unsteady form, Eqs. (4.17) and (4.19), does not yield better correlations between the 
computed and experimental results. Therefore, the steady form, Eqs. (4.20), was used 
also for unsteady flow cases. Further studies are needed to clarify whether the explicit 
time integration of Eqs. (4.15) is consistent with the implicit Navier-Stokes and Full- 
Potential solvers used in the present work. 

With the resulting values of ci to cs> the changes in flow properties at the 
interface are computed using the inverse of Eq. (4.14): 
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4.3.2. Interface Boundary Conditions for the Full-Potential Solver 


From the last two equations in Eq. (4.14): 


c --2 


(4.22) 


Using 


p = ±(a^=>5p = ^ [ Sa 


(4.23) 


yields: 


CA 


pa 

2^[A6 


5W + 


8a 


Cs = -&Lsw + 
24 m y-i 


Pa 

y-1 
pa Sa 


(4.24) 


By multiplying Eq. (4.24) by 2/ pa the Riemann invariants are recovered, i.e., 
the characteristic invariants may be expressed as: 


SW . 2 

Ri = — ;= + 

VS r- 1 

51V 2 . 

Rl — r—~ "h . 

Vd6 r-i 


(4.25) 
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The first characteristic corresponds to the eigenvalue A 4 = W-a^pU, and the 
second to As = W- + a^Ae . The Riemann invariants are computed according to the 
signs of the corresponding eigenvalues, for example if A 4 > 0 then R x is computed 
using information from the inner zone, otherwise it is computed using information from 
the outer zone. With the resulting values of R x and R 2 , the changes in flow properties 
at the interface are computed using the inverse of Eq. (4.25). 

As shown in Chapter V , the above procedure has been successful in suppressing 
the oscillatory behavior observed in the computations with the original boundary 
conditions. Although the procedure is strictly valid only for steady flows, it was used 
also for unsteady flows with results similar to those obtained by full Navier-Stokes 
computations, as discussed in Chapter V. 
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CHAPTER V 


RESULTS AND DISCUSSION 


In the present Chapter, applications of the hybrid Navier-Stokes/Full-Potential 
Method are discussed. First, the application of the method to a rectangular wing in 
steady flow is presented. Next, the application to an F-5 wing in steady and unsteady 
flow for both subsonic and transonic flow is presented. The method was also applied to 
a rectangular wing with a simulated glaze ice accretion in subsonic flow, and the 
corresponding results are presented in Appendix A for completeness. 

5.1. Rectangular Wing Study 

5.1.1. Configuration 

The hybrid Navier-Stokes/Full-Potential Method has been applied to a 
rectangular wing of aspect ratio 2.5. The airfoil section was NACA 0012. This 
configuration has been experimentally studied by Bragg et al. 151 " 153 as part of their 
iced wing studies. The surface pressures were measured at five spanwise stations: 17%, 
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34%, 50%, 66% and 85% of the wing semi-span. The wind tunnel had a 0.85m by 
1.22m rectangular cross section. The model semi-span was 0.95m. 



Fig. 5.1: Computational Grid for the Rectangular NACA 0012 Wing. 


The computational grid used in the present study, illustrated in Fig. 5.1, was an 
algebraic C-grid with 141 x 19 x 41 grid points, with 121 points over the airfoil surface 
at each span wise station. 14 spanwise stations were used along the wing, with 5 stations 
extending beyond the tip. The Navier-Stokes and Full-Potential solvers were interfaced 
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at KMATCH— 20, so that about half of the number of points are located in each zone. 
The Mach number was 0.12, the Reynolds number 1.5 million, and the angles of attack 
were 4° and 8°. The CPU time for this configuration was about 13.6 sec. per iteration 
on a Hewlett-Packard Apollo 700 workstation. This time was about 60% of the time 
required for a full N avier-Stokes analysis. The CPU time for the same configuration on 
Georgia Tech’s Cray Y-MP/E was about 4.5 sec. per iteration. This time was about 68% 
of the time required for a full Navier-Stokes analysis. It was observed that the Navier- 
Stokes module vectorized better than the Full-Potential module, therefore the 
computational savings on a vector processor were lower. 

5.1.2. Surface Pressure Distributions 

The pressure coefficient distributions for angles of attack of 4° and 8° are 
shown in Figs. 5.2 and 5.3, respectively. The small spurious peaks near the trailing 
edge at all the spanwise stations are due to inadequate grid resolution in that region, 
where the chordwise spacing was about 4%, and may be improved with a clustered 
grid. Outside the trailing edge region, a very good agreement may be observed for the 
case (X— 4 , with the possible exception of station 85%. The discrepancies at that station 
may be due to the tapering of the grid outboard of the wing tip, and better results could 
be obtained with a finer grid near the tip. In Fig. 5.3, corresponding to the case a=8°, 
the results obtained from full Navier-Stokes computations are shown for comparison. 
For this case, an underprediction of the suction peak is noticeable. It can be observed 
that the undeiprediction of the suction peak for the case a=8° also occurs for the full 
Navier-Stokes computations. It should be noted that the experimental results used here 
were obtained with a small clearance between the upper and lower wind-tunnel walls 
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and the wing surface. Consequently, wind-tunnel wall interference may partially 
account for the underpiediction of the suction peaks. 

Further insight can be obtained by analyzing the lift coefficient distribution 
along the span, shown in Fig. 5.4. The underprediction of lift can be clearly noticed in 
this figure. Overall, it is observed that our current results correlate well with 
experiments and are consistent with those obtained by a full Navier-Stokes code 8 " 10 , 
while consuming only about 60% of the CPU time. 
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Fig. 5.2: Surface Pressure Distributions for the Rectangular Wing at a=4° 
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Fig. 5.2: Surface Pressure Distributions for the Rectangular Wing at a=4° (continued) 
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Fig. 5.2: Surface Pressure Distributions for the Rectangular Wing at a=4° (concluded) 



Fig. 5.3: Surface Pressure Distributions for the Rectangular Wing at a=8° 
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Fig. 5.4: Lift Distribution for the Rectangular Wing at a=4 and 8 C 


5.1.3. Convergence Histones 

The convergence histories for the original hybrid scheme showed a highly 
oscillatory behavior, as shown in Fig. 5.5 for the rectangular wing at a=8°. These 
oscillations indicated false reflections from the Navier-Stokes/Full-Potential interface, 
which affected the solution at early time levels, and required a large number of 
iterations, about 4000 to 6000, while the Navier-Stokes code was able to achieve 
satisfactory convergence in about 2000 iterations. 
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A significant improvement was obtained by using the characteristic-based 
interface boundary conditions described in Chapter IV. The convergence history of the 
maximum residual in the Navier-Stokes zone for the rectangular wing at a=4° is shown 
in Fig. 5.6. The convergence history of the full Navier-Stokes computation is also 
shown for comparison. It can be observed that the hybrid method attains convergence 
characteristics similar to the full Navier-Stokes method when proper interface boundary 
conditions are used. Consequently the hybrid method was able to fully realize the 
computational savings of about 40% when the improved boundary conditions were 
used. 



Iteration 


Fig. 5.5: Residual Histories in Potential and Navier-Stokes Zones 
(Rectangular Wing at a=8°). 
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Max Residual 



Fig. 5.6: Residual History in Navier-Stokes Zone with Original and Improved 
Interface Boundary Conditions (Rectangular Wing at a=4°). 



5.2. F-5 Wing Studies 


In order to investigate the applicability of the hybrid Navier-Stokes/Full- 
Potential method to unsteady compressible viscous flows, an unsteady problem 
involving the F-5 wing in pitch oscillations was studied, so that the effect of the 
interface boundary conditions on the accuracy of the time-dependent results could be 
verified. The investigation presented here is aimed at the transonic range, which is a 
very rigorous test for the present method, due to the development of shock waves 
which cross the Navier-Stokes/Full-Potential interface. The interface boundary 
conditions are therefore required to propagate significant disturbances. In the unsteady 
flow simulations, these disturbances have to be propagated in a time-accurate fashion, 
which presents an even more rigorous test. However, unsteady transonic flow 
simulations present also a potential for significant computational savings by using the 
present method, because of the numerous computations that are needed in a typical 
parametric investigation. 


5.1.1. Configuration 

The F-5 wing is chosen because a wealth of steady and unsteady data, as well as 
detailed geometric description, are readily available. It also represents a challenging 
configuration, since features such as taper and a thin, drooped leading edge are present 
The experimental lay-out used by Tijdeman et al 43 * 44 is illustrated in Fig. 5.7. The 
investigators measured steady and time-dependent pressures at eight spanwise stations 
indicated in Fig. 5.7 and listed in Table 5.1. Note that, for all cases presented here, no 
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experimental data are used for the first point in the upper surface of the wing at 
spanwise stations 3 and 5, because the measuring probes at these points were defective. 



Fig. 5.7: F-5 Wing Experimental Lay-Out 
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Table 5.1: Spanwise Stations Where Experimental 
Data is Available for the F-5 Wing 


— 

Station No. 

■SI 

1 

18.1 

2 

35.2 

3 

51.2 

4 

64.1 

5 

72.1 

6 

81.7 

7 

87.5 

8 

97.7 


The computational grid used in the present study, illustrated in Fig. 5.8, was an 
algebraic C-grid with 141 x 19 x 41 grid points, with 121 points over the wing surface 
at each spanwise station. 14 spanwise stations were used along the wing, with 5 stations 
extending beyond the tip. The Navier-Stokes and Full-Potential solvers were interfaced 
at KMATCH=20, so that about half of the number of points are located in each zone. 
The computations presented here were performed in NASA Lewis Research Center 
Cray-Y/MP. The CPU time for this configuration was about 0.95 sec per time step, 
about 73% of the CPU time needed for a full Navier-Stokes computation. 
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Fig. 5.8: Computational Grid for the F-5 Wing 


5.2.2. Steady Flow Simulations 

For steady flow simulations, the Mach numbers used in the present study were 
0.6, 0.8, 0.9 and 0.95. The Reynolds number based on the root chord and fxee-stream 
speed = pJJ~,CrlV'~ was 11 million for Moo=0.8. The angle of attack for all 
computations presented here was 0°. 
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The initial test case was Moo=0.6. The steady pressure coefficient distributions 
at spanwise stations 35.2%, 72.1% and 97.7% are shown in Figs. 5.9(a), (b), and (c), 
respectively. The numerical results are in very good agreement with the experimental 
data. In particular the leading edge suction peaks seem well predicted, although 
insufficient experimental data are available in that region. Also the precise location 
where the lower and upper surface pressure coefficients match seems well predicted. 

The next test case was Moo=0.8. The steady pressure coefficient distributions at 
spanwise stations 35.2%, 72.1% and 97.7% are shown in Figs. 5.10(a), (b), and (c), 
respectively. Again the main features of the pressure distributions are well predicted. 

A more demanding test case was Moo=0.9, for which the steady pressure 
coefficient distributions at spanwise stations 35.2%, 51.2%, 64.1%, 72.1%, 87.5% and 
97.7% are shown in Figs. 5.11(a) through (f), respectively. Again the suction peaks and 
the location where upper and lower surface pressure coefficient match seem well 
predicted. A slight underprediction of the suction values between 10% and 50% of the 
chord on the upper surface at the inboard stations is apparent. But the main feature of 
this configuration is the onset of a weak shock on the upper surface, near the tip. This 
shock may be observed in the experimental data in Fig. 5.11(f), but it appears that the 
weak shock was smeared in the numerical solution and cannot be clearly identified. By 
checking for supersonic points in the full-potential region, it was observed that the 
sonic line did cross the Navier-Stokes/Full-Potential interface for this configuration. 
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The last steady test case was Moo=0.95, for which the steady pressure coefficient 
distributions at spanwise stations 35.2%, 51.2%, 64.1%, 72.1%, 87.5% and 97.7% are 
shown in Figs. 5.12(a) through (f), respectively. Here the dominating feature is the 
shock that forms over most of the wing, on both upper and lower surfaces. The upper 
surface shock is stronger and aft of the lower surface shock. These features were well 
predicted by the current method, although some slight underprediction of the mid-chord 
upper surface suction values is noticed. The suction peaks and location of matching 
upper and lower surface pressures are again well predicted, except at the station 97.7%, 
where the experimental data indicate a lower suction peak. This test case was 
demanding in the sense that the shock crosses the Navier-Stokes/Full-Potential 
interface, and the results presented here indicate that the hybrid method is able to 
predict adequately both shock location and strength even when the discontinuities due 
to the shock are propagated through the Navier/Stokes/Full-Potential interface. Further 
evidence to support this conclusion is presented in Fig. 5.13, where the density 
contours at station 81.7% of span, Moo = 0.95, are shown. In this figure, the Navier- 
Stokes/Full-Potential interface is drawn to facilitate the analysis. It can be seen that the 
contours smoothly cross the interface, and in particular the shock is well captured 
across the interface. 

Overall the results presented here show that the hybrid method can be 
successfully applied to steady transonic flows, even when the shock crosses the Navier- 
Stokes/Full-Potential interface. The discrepancies observed near the tip seem to be due 
to inherent inaccuracies associated with the tapering of the computational grid outside 
the tip, and are present also in full Navier-Stokes computations, and consequently not 
associated with deficiencies in the hybrid method. Furthermore, the differences 
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between the computed and measured data are of the same order of, or lower than, the 
more costly Navier-Stokes solutions presented by other researchers 30 ’ 88 * 154 . Numerical 
experimentation also indicated that for the transonic cases presented here, the current 
hybrid code with characteristic-based interface boundary conditions achieved higher 
convergence rates than the original hybrid code and even the full Navier-Stokes 
computations, probably due to the use of a higher time step in the FPE solver.Thus, it 
may be concluded that the present hybrid method gives results that are comparable to 
the more exact approaches, at less than 60% of the cost of Navier-Stokes simulations. 
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(c) 97.7% Span 

Fig. 5.9: Steady Surface Pressure Distributions, Moo = 0.6 (concluded) 
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Fig. 5.10: Steady Surface Pressure Distributions, Moo = 0.8 
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(c) 97.7% Span 

Fig. 5.10: Steady Surface Pressure Distributions, Moo = 0.8 (concluded) 
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(d) 72.1% Span 

Fig. 5.11: Steady Surface Pressure Distributions, ML, = 0.9 (continued) 





(f) 97.7% Span 

Fig. 5.1 1: Steady Surface Pressure Distributions, Me = 0.9 (concluded) 
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Fig. 5.12: Steady Surface Pressure Distributions, = 0.95 
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5.23. Unsteady Flow Simulations 


Numerous unsteady flow simulations are typically needed to predict aeroelastic 
characteristics. These predictions are especially important in transonic flow, where 
aeroelastic instabilities are more likely to occur. In the transonic regime, the flow is 
inherendy nonlinear, which has prevented the use of simpler methods, such as panel 
methods and vortex lattice methods. In particular, the presence of a supersonic region 
embedded in a subsonic region causes downstream disturbances to be propagated 
upstream with a considerable time lag, which results in significant out-of-phase forces. 
Full-potential and Euler methods have been applied to such flows with some success, 
but viscous effects can alter the location and strength of shocks. In an unsteady flow 
situation, both the location and strength of the shocks can change rapidly and generate 
significant unsteady forces, which have to be known for aeroelastic analysis. Navier- 
Stokes analyses are therefore a natural choice, but require drastically high 
computational resources. The present hybrid Navier-Stokes/Full-Potential method 
presents the advantage of fully capturing the viscous and nonlinear effects, while 
incurring a significantly lower computational cost. However, unsteady transonic flows 
are also very challenging to the present method, due to the presence of strong 
disturbances generated by unsteady shock motion, which need to be propagated 
through the Navier-Stokes/Full-Potential interface in a time-accurate fashion. The 
purpose of the investigation presented in this Section is to validate the present method 
for these flow conditions. 
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For the unsteady flow simulations presented here, the Mach numbers were 0.6, 

0.9 and 0.95. The wing was in pitching oscillations about half-chord, and around 
oco=0°> i.e.: 


a(r) = ao + Aa sin(xr) 


(5.1) 


where Aoc is the amplitude of oscillation, given in Table 5.2 for the various test cases. 
Also, T — a~tjc is the nondimensional time, and K = (ocja - is the reduced frequency, 
with c the reference length. Note that co = 2 tzF , where F is expressed in Hertz. Note 
also that the definition of reduced frequency here differs from that of Tijdeman et al. 43 , 
for consistency with the nondimensionalization used in the present work. The definition 
of reduced frequency here may be related to the definition in the cited reference as: 


K 


(PC 

doo 


2itFc 

do a 


nFcr Uoo 2c 

FJ oo flw Cr 



Cr 


(5.2) 


where Ki is the reduced frequency as defined by Tijdeman et al. and Cr is the root 
chord, equal to 0.6396 meters. The reference length used in the present study was 1 

meter. The reduced frequencies corresponding to the test cases used here are listed in 
Table 5.2. 
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Table 5.2: Test Cases for Unsteady Flow 


Test Case 


ess 


mm 

mmmmm 


0.6 

20 

0.106 

0.199 

0.373 


0.6 

40 

0.115 

0.399 

0.749 


0.9 

40 

0.111 

0.275 

0.774 

i i 

0.95 

40 

0.222 

0.264 | 

0.785 


Under these pitching oscillations, the F-5 wing deforms aeroelastically. During 
the investigation reported by Tijdeman et al.^3,44^ the wing vibration mode was 
measured for the various test runs using eight accelerometers. These measurements 
were used to obtain an approximate analytical expression for the vertical wing 
displacement at various points, assuming no chordwise deformation (rigid rotation) and 
parabolic spanwise deformation: 

w ( x *y) ~ «oo + ao\X + aioy + an xy + a^y 2 + +a 2 ixy 2 (5.3) 

where the coefficients a i} are tabulated in Ref. 44. This approximation to the elastic 
deformation allows a consistent representation by a rigid rotation about the node 
corresponding to each spanwise station. The nodal lines corresponding to the cases 
presented here are illustrated in Fig. 5.14, from Ref. 44. 
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(c) Moo=0.9, F=40Hz (d) IVU=0.95, F=40Hz 

Fig. 5.14: Nodal Lines forF-5 Pitching Oscillations 


For a given spanwise station yj, the node is located at x^j and the motion at that 

spanwise station is represented by a local rigid pitching rotation about this node, given 
by: 


(Xj(T) = 9jCc(T) ( 5 . 4 ) 

where a(r) is given by (5.1) and Qj represents the mode shape, obtained from (5.3). 
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For each unsteady flow simulation presented here, the job is restarted from a 
previously converged steady flow solution at the specified Mach number and average 
angle of attack, <Xo=0° in all cases considered here. It should be noted that the response 
of the flow field to the oscillatory motion will involve a transient and a steady-state 
response. This is illustrated in Fig. 5.15, where the time history of the lift coefficient at 
the spanwise station 64.1% for the case Moo=0.95, F=40 Hz, is shown. The results that 
will be presented subsequently involve only the steady-state response. For that purpose, 
enough iterations — typically about one cycle — were ran before the computation of 
the unsteady pressure coefficients was started. For the case Mo<>=0.9, F=40 Hz, an 
additional cycle was ran and the unsteady pressure coefficients were recomputed, 
yielding almost the same results with no perceptible differences from the previous 
cycle. This indicated that one cycle was indeed enough to eliminate the transient. 
Another concern was to verify that the response was dominated by the first harmonic. 
For that purpose, the second harmonic was also computed and it was observed that its 
values were negligible, except at locations where the first harmonics presented peaks, 
but at these locations the second harmonics were still much smaller than the first 
harmonics. 
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Fig. 5.15: Time History of the Lift Coefficient at 64.1% Span, Moo=0.95, F=40 Hz. 

The unsteady pressure coefficients were computed as follows: The actual 
oscillations of the wing are given by (5.2) with (Xo=0°. Let the complex representation 
of the motion be: 


a c (r) = Aa e iKC 


(5.5) 


Now let the complex representation of the steady-state response be: 


Cp(t) = [Re(Cj>j) + /lm(0/)] e‘" 


(5.6) 


88 




Since the actual motion is given by the imaginary part of (5.5), the actual 

response is given by the imaginary part of (5.6) and can be related to the complex 
representation by: 

Cr(t)= [Re(Cpi) + i Im(C/.j)] e'" 

= [Re(C P ,)cos(xr) - Im(Cm)sin(Kr)] + i[Re(Cp,)sin()cr) + Im(C/>,)cos(xr)] (5.7) 

CM ’ 

Then Re(c Pl ) and Im(Cf,) can be obtained in terms of the actual pressure 
coefficient response by: 

jg ti+2u/ic 

Re(o t ) = — jCp( t) sin( /rr)dT 

n tt 

]£ ti+2xfK (5.8) 

Im(C/>i) = — j Cp(z) cos(m)dT 
n Tl 


where ii is chosen so that the transient is not included in the computations. It should be 

noted that the experimental unsteady pressure coefficients 43 * 44 were normalized with 
respect to the amplitude of oscillations, 2Aa. The actual computation of Re(C/>,) and 
Im(Cpi) is finally performed by the discretization of Eq. (5.8): 

_ / v kAt m ' +m r 

Re(Cw)= 2^if- sin(wnAT ) 

T s v kAt nti+mr (5.9) 

MCpi)= 2^; J^ C0S{KmA ^ 
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where m denotes the time step number, rm = Ti/At, and m T = 2n/K^i. ;. Note that the 

normalization with respect to the amplitude of oscillations has been included in Eq. 
(5.9). 


Using the above described procedure, the unsteady pressure coefficients were 
computed for the configurations listed in Table 5.2. The first test case corresponds to 
Mach number 0.6, frequency of 20 Hz and amplitude of oscillation 0.106 degrees. The 
nondimensional time step AT=a„Ar/c was 0.005, so that 3369 time steps per cycle 
were needed. The unsteady pressure coefficient distributions at spanwise stations 
35.2% and 72.1% are shown in Fig. 5.16. The in-phase (real) component is in good 
agreement with the experiment, although few test data are available near the leading 
edge to confirm the level of unsteady pressure peaks. In particular, at 72.1% span the 
first experimental point in the upper surface is missing and the first experimental point 
in the lower surface seems faulty, as a peak is expected near the leading edge. The out- 
of-phase (imaginary) component leading edge peaks appear well predicted, and so are 
the distributions from about 40% chord towards the trailing edge, although the location 
where the lower and upper surface imaginary components match is predicted aft of the 
experimental location. 

The next test case corresponds to Mach number 0.6, frequency of 40 Hz and 
amplitude of oscillation 0.115 degrees. The nondimensional time step was 0.005, so 
that 1678 time steps per cycle were needed. The unsteady pressure coefficient 
distributions at spanwise stations 35.2% and 72.1% are shown in Fig. 5.17. Again the 
in-phase (real) component appears to be in good agreement with the experiment, 
showing little change with respect to the previous test case (F=20Hz). At this higher 
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frequency, the out-of-phase (imaginary) component is overall higher than for the lower 
frequency, in both the experimental and computed results. This imaginary component 
appears well predicted for the locations aft of 30% chord, but the location where the 
lower and upper surface imaginary components cross is predicted aft of the 
experimental location. The leading edge peaks also seem overpredicted, although very 
few experimental points are available there. 

The next test case was Moo=0.9, frequency 40 Hz and amplitude of oscillation 
0.111 degrees. The nondimensional time step was 0.005, so that 1624 time steps per 
cycle were needed. The real and imaginary parts of the unsteady pressure coefficient 
distributions at spanwise stations 35.2%, 51.2%, 64.1%, 72.1%, 87.5% and 97.7% are 
shown in Figs. 5.18(a) through (1). At all spanwise stations strong leading edge peaks 
are present in the lower surface both in the in-phase and out-of phase component and 
they seem to be well predicted by the present method. It may be recalled that in the 
steady case, discussed in Section 5.2.2 above, the main feature of this configuration 
was the onset of a weak shock on the upper surface, near the tip, which could be 
observed in the experimental data in Fig. 5.1 1(f). In the unsteady case, strong variations 
in the pressure coefficients are noticed about 50% of the chord at all spanwise stations. 
These indicate that the oscillatory motion causes a shock wave to form over most of the 
wing and move back and forth, generating the strong variations observed here. 


On the lower surface, the real part of the unsteady pressure coefficient shows 
relatively small variations, as seen in Figs. 5.18(a), (c), (e), (g), (i) and (k), while the 
imaginary part shows stronger variations. This indicates that the shock wave on the 
lower surface is predominantly out of phase with respect to the motion. The present 
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method predicts this phenomenon for inboard stations, up to about half span, as seen in 
Figs. 5.18(b) and (d), although about 15% of the chord aft of the experimental data and 
somewhat smeared. At outer stations, the unsteady shock on the lower surface appears 
to smeared and even disappears in the computed results. 

On the upper surface, both the real and imaginary components show significant 
variations indicating a stronger shock moving, with a significant in-phase, but higher 
out-of-phase component. This phenomenon is better predicted in all spanwise stations, 
especially the chordwise location, around mid-chord. The imaginary part in Figs. 
5.18(b), (d), (f) and (h) seems to indicate some underprediction of strength of the out- 
of-phase component of the shock, which may be due to smearing. 

The results discussed above for the case Moo=0.9, F=40 Hz are quite similar to 
full Navier-Stokes computations performed with a different version of the Navier- 
Stokes module used here 30 . The same test case was also investigated by Obayashi et 
al. 154 , who used a streamwise upwind algorithm. In that paper, only the upper surface 
results were presented and a comparison was made between the streamwise upwind 
algorithm and a central-difference algorithm, with a more favorable correlation for the 
former. Since the present method is modular, an upwind Navier-Stokes module could 
potentially improve the unsteady shock prediction while maintaining the computational 
savings obtained by the hybrid Navier-Stokes/Full-Potential method. It should be noted 
that this case was a rigorous test for the hybrid method, since the unsteady shock is 
moving back and forth and crossing the Navier-Stokes/Full-Potential interface. As the 
shock is not completely aligned with the grid, strong oblique disturbances are 
propagated through the interface. The results presented here show that the interface 
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boundary conditions are adequate to transmit these disturbances in a time-accurate 
fashion. 

The last test case was Moo=0.95, frequency 40 Hz and amplitude of oscillation 
0.222 degrees. The nondimensional time step was 0.005, so that 1601 time steps per 
cycle were needed. An additional computation to assess the effect of the time step was 
made with a nondimensional time step of 0.002, for which 4003 time steps per cycle 
were needed. The real and imaginary parts of the unsteady pressure coefficient 
distributions at spanwise stations 35.2%, 51.2%, 64.1%, 72.1%, 87.5% and 97.7% are 
shown in Figs. 5.19(a) through (1). As occurred in the previous test case, at all spanwise 
stations strong leading edge peaks are present in the lower surface both in the in-phase 
and out-of phase component and they seem to be well predicted by the present method. 
The steady flow results for Moo=0.95, shown in Fig. 5.12, indicate a strong shock on 
both upper and lower surfaces around 80% of the chord. The experimental data for the 
unsteady case, seen in Fig. 5.19, show significant peaks around this chordwise location, 
mostly in the real (in-phase) component, but also in the imaginary (out-of-phase) 
component. These peaks are very localized, which indicates that they result more from 
shock strength variations than shock movement The numerical results presented in Fig. 
5.19 show that the present method was unable to correctly predict the peak in the real 
part, but predicted the peak in the imaginary part. The computations with a smaller time 
step show some improvement in the real part, which indicate that the time step might 
have to be further reduced to yield a better correlation. Further reductions in time step 
were not attempted because of the large CPU resources that would be needed. This 
improvement with time step also indicates that the deficiency is not due to the hybrid 
method. It should also be noted that the current coarse grid presents some smearing in 


93 



the shock, therefore small changes in the shock strength are not likely to be well 
captured, even with a smaller time step. 

Except for the above discussed discrepancy, the unsteady pressure coefficient 
distribution is well predicted. It should be noted that this is a very rigorous test for the 
present method, due to the strong shock crossing the Navier-Stokes interface, as seen in 
Fig. 5.13. The results presented here indicate that the discrepancies observed in this test 
case are inherent to the Navier-Stokes module, and can probably be overcome by using 
an upwind Navier-Stokes module capable of capturing sharper shocks. 

Overall the unsteady pressure coefficient distributions correlate well with 
experimental data and are similar to those obtained with equivalent full Navier-Stokes 
computations, with a fraction of the computational cost. The savings in CPU time were 
found to depend on the vector capability of the CPU, ranging from 27% on the Cray 
Y/MP-L up to 40% on a HP Apollo 700 workstation. It is believed that effective clock 
times could be reduced even more with respect to full Navier-Stokes computations on 
distributed processing machines, since the entire Full-Potential module could be solved 
in parallel without data exchange with the Navier-Stokes solver, for a given iteration. 
This possibility, however, is yet to be investigated. 
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(b) Imaginary Part, 35.2% Span 

Fig. 5.16: Unsteady Surface Pressure Distributions, Me*, = 0.6, f=20 Hz 
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(d) Imaginary Part, 72.1% Span 

Fig. 5.16: Unsteady Surface Pressure Distributions, Moo = 0.6, f=20 Hz (concluded) 
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(b) Imaginary Part, 35.2% Span 

Fig. 5.17: Unsteady Surface Pressure Distributions, M«, = 0.6, f=40 Hz 


97 




(d) Imaginary Part, 72.1% Span 

Fig. 5.17: Unsteady Surface Pressure Distributions, Moo = 0.6, f=40 Hz (concluded) 
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(d) Imaginary Part, 51.2% Span 

Fig. 5.18: Unsteady Surface Pressure Distributions, Moo = 0.9, f=40 Hz (continued) 
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(f) Imaginary Part, 64.1% Span 

Fig. 5.18: Unsteady Surface Pressure Distributions, M«> = 0.9, f=40 Hz (continued) 








(h) Imaginary Part, 72.1% Span 

Fig. 5.18: Unsteady Surface Pressure Distributions, Moo = 0.9, f=40 Hz (continued) 
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(j) Imaginary Part, 87.5% Span 

Fig. 5.18: Unsteady Surface Pressure Distributions, Moo = 0.9, f=40 Hz (continued) 


103 







(1) Imaginary Part, 97.7% Span 

Fig. 5.18: Unsteady Surface Pressure Distributions, Moo = 0.9, f=40 Hz (concluded) 
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(a) Real Part, 35.2% Span 



(b) Imaginary Part, 35.2% Span 

Fig. 5.19: Unsteady Surface Pressure Distributions, M* = 0.95, f=40 Hz 
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(d) Imaginary Part, 51.2% Span 

Fig. 5.19: Unsteady Surface Pressure Distributions, Moo = 0.95, f=40 Hz (continued) 
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(f) Imaginary Part, 64.1% Span 

Fig. 5.19: Unsteady Surface Pressure Distributions, J4 = 0.95, f=40 Hz (continued) 
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(h) Imaginary Part, 72.1% Span 

Fig. 5.19: Unsteady Surface Pressure Distributions, = 0.95, f=40 Hz (continued) 
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(i) Real Part, 87.5% Span 



(j) Imaginary Part, 87.5% Span 

Fig. 5.19: Unsteady Surface Pressure Distributions, M*, = 0.95, f=40 Hz (continued) 
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Fig. 5.19: Unsteady Surface Pressure Distributions, Moo = 0.95, f=40 Hz (concluded) 
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CHAPTER VI 


CONCLUSIONS AND RECOMMENDATIONS 


6.1. Conclusions 

An improved hybrid Navier-Stokes/Full-Potential method has been developed 
This method has been applied to subsonic and transonic steady and unsteady flow 
cases. The following conclusions can be drawn from the present investigation: 

• The present method is an economical way of studying steady and unsteady 
flows. The accuracy of Navier-Stokes computations can be retained with savings in 
CPU times of more than 40% in some cases. 

• When coupled analysis are used, special attention must be paid to interface 
boundary conditions. Improper boundary conditions allow false reflections, which can 
slow down convergence to steady-state or lead to loss of temporal accuracy. It was 
found that the implementation of characteristic-based interface boundary conditions 
developed in the present work can adequately treat the interface, and allow signals to 
pass to and from one zone to another. This is especially important when shocks are 
present, because strong oblique disturbances must be transmitted through the interface. 


Ill 


• The interface boundary conditions developed in the present work are strictly 
valid only for steady flows. However, it has been found that they could be applied to 
unsteady flow cases with results similar to those obtained by full Navier-Stokes 
computations. 

• For unsteady transonic flows, Navier-Stokes computations are clearly needed, 
because the strength and location of shocks are a major factor in determining unsteady 
loads. Full-Potential methods predict shocks aft of their actual location and overpredict 
their strength. 

• The savings in CPU time were found to depend on the vector capability of the 
CPU, ranging from 27% on the Cray Y/MP-L up to 40% on a HP Apollo 700 
workstation. For steady flow cases, the computational savings were slightly higher, 
because the hybrid method presented convergence rates higher than the full Navier- 
Stokes method. 


6.2. Recommendations 

• The present hybrid Navier-Stokes/Full Potential method has proven to be an 
effective way to maintain the accuracy of the Navier-Stokes simulations with 
substantial reductions in computational cost The savings in CPU time were found to 
depend on the vector capability of the CPU, ranging from 27% on the Cray Y/MP-L up 
to 40% on a HP Apollo 700 workstation. For this reason, it is recommended that further 
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improvements in vector processing machines be pursued by further vectorization of the 
FPE solver. 

• Effective clock times could be reduced even more with respect to full Navier- 
Stokes computations on distributed processing machines, since the entire Full-Potential 
module could be solved in parallel without data exchange with the Navier-Stokes 
solver, for a given iteration. It is recommended that the implementation of the present 
hybrid method on such machines be investigated. 

• It is recommended that further studies be conducted to ascertain the suitability 
of the interface boundary conditions developed in the present work to unsteady flows. 
These studies should include an investigation on the consistency of applying explicit 
time integration at the viscous/inviscid boundary, while using implicit solvers in the 
inner and outer regions. 

• The shock capturing capability of the present method seems to be limited by 
the computational grid, in the inner region. The test case of the F-5 wing in pitching 
oscillations at Moo=0.95, with a strong shock moving over the aft part of the wing, 
presents an opportunity for grid sensitivity studies. It is recommended that these 
sensitivity studies are performed, to ascertain whether the inability of the present 
method to predict the peak in the in-phase component of the unsteady pressure 
coefficients is due solely to the lack of adequate grid refinement in the aft part of the 
wing. 


113 



• Since the present method is modular, numerical experiments with different 
solvers could be undertaken. It is recommended that the ensemble FPE solver + 
interface boundary conditions be constructed as a ‘plug-in* module to be applied to 
existing Navier-Stokes solvers as a means for quickly reducing the computational cost 
of existing Navier-Stokes methods. 

• It is recommended that the present methodology be studied in connection with 
unstructured grid - based Navier-Stokes solvers. These methods are computationally 
intensive and are likely to benefit substantially from the computational savings allowed 
by the present technique. 

In closing, it is hoped that the present hybrid technique, which combines the 
accuracy of Navier-Stokes methods in the viscous regions with the economy of 
potential flow methods in inviscid regions, will be used as a stepping stone for more 
ambitious efforts involving aeroelastic and unsteady aerodynamic analysis of complete 
aircraft configurations. 
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APPENDIX A 


APPLICATION OF THE HYBRID NAVIER- 
ST OKES/FULL-POTENTIAL METHOD TO AN 

ICED WING 


In this Appendix, the application of the hybrid Navier-Stokes/Full-Potential 
Method to a rectangular wing with a simulated glaze ice accretion in steady flow is 
presented. 


A.l. Configuration 

The hybrid Navier-Stokes/Full-Potential Method has been applied to an unswept 
wing of aspect ratio 2.5 with a simulated glaze ice accretion as shown in Fig. A.l. This 
configuration has been experimentally studied by Bragg et al.l 5 ! -153 . The results for 
the same configuration without the simulated ice accretion were presented in Chapter 
V. The surface pressures were measured at five spanwise stations: 17%, 34%, 50%, 
66% and 85% of the wing semi-span. 
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Fig. A.l: Simulated Glaze Ice Accretion 


The computational grid used in the present study was an algebraic C-grid with 
141 x 19 x 41 grid points, with 111 points over the airfoil (and ice) surface at each 
spanwise station. 14 spanwise stations were used along the wing, with 5 stations 
extending beyond the tip. The computational grid is illustrated in Fig. A.2. The Navier- 
Stokes and Full-Potential solvers were interfaced at KMATCH=20, so that about half 
of the number of points are located in each zone. The Mach number was 0.12, the 
Reynolds number 1.5 million, and the angles of attack were 4° and 8°. The CPU time 
for this configuration was about 24 sec. per iteration on a Hewlett-Packard Apollo 700 
workstation. The CPU time for the same configuration on Georgia Tech's Cray Y-MP/E 
was about 9.6 sec. per iteration. These times are about 60% of the times required for a 
full Navier- Stokes analysis. 


A*2« Surface Pressure Distributions 

The computed surface pressure distributions at the computational span stations 
were linearly interpolated to the span stations where experimental data is available to 
allow direct one-to-one comparisons^ 1-153 # The resulting pressure coefficient 
distributions for the angles of attack 4° and 8° angle of attack are shown in Fig. A.3 and 
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A.4, respectively. The spurious peaks near the trailing edge at all the spanwise stations 
are due to inadequate grid resolution in that region. A good agreement may be observed 
between the numerical and experimental results, except for the spurious peaks near the 
trailing edge as mentioned above and for the strong peaks and oscillations near the 
separated region. It should be noted that essentially the same results are obtained with 
Navier-Stokes computations. For the 4° angle of attack configuration, the length of the 
separation bubble appears very well predicted, although the exact chordwise location 
where the it starts appears displaced, except near the tip. Due to the sharp pressure 
variations just before the separation bubble, this slight displacement causes the pressure 
level along the separated region to be somewhat underpredicted. For the 8° angle of 
attack configuration, the displacement of chordwise location where the separation 
bubble starts is observed only inboard, but the length of the separation bubble is not 
predicted as well as for the 4° configuration. 

It should be noted that the experimental results used here were obtained with a 
small clearance between the upper and lower wind-tunnel walls and the wing surface. 
Consequently, wind-tunnel wall interference may partially account for the 
discrepancies noted above. 

Further insight can be obtained by analyzing the lift coefficient distribution along 
the span, shown in Fig. A.5. Interestingly, the underprediction of lift due to the lower 
suction peak, observed for the clean wing configuration, is not observed for the 4° 
angle of attack iced configuration. This may be attributed to the sharp variation of 
pressure coefficient with formation of the separation bubble near the leading edge. 
Since the suction peak occurs in a much narrower region than in the non-iced 
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configuration, the actual peak value becomes less important for the section lift The 
underprediction of lift at 8° angle of attack iced configuration at inboard stations, in 
turn, appears to be due to the inaccurate prediction of the start of the separation bubble, 
associated with its relatively longer length for this configuration. 

Overall, it is observed that our current results correlate well with experiments and 
are consistent with those obtained by a full Navier-Stokes code 8 - 10 , while consuming 
only about 60% of the CPU time. 



Fig. A.2: Computational Grid for Rectangular Wing with Simulated Glaze Ice 

Accretion 
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(b) 34% Span 

Fig. A.3: Surface Pressure Distributions for Iced Wing at a=4° 















Fig. A.4: Surface Pressure Distributions for Iced Wing at a=8° (continued) 
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(e) 85% Span 

Fig. A.4: Surface Pressure Distributions for Iced Wing at a=8° (concluded) 







Fig. A.5: Lift Distribution for Iced Wing at a =4 and 8° 
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